## Code to replicate Table 3 of Abrajano, Elmendorf, and Quinn. "Measuring
## Perceived Skin Color: Spillover Effects and Likert-type Scales". JOP.


library(clubSandwich)


sink(file="AbrajanoElmendorfQuinn-Table3-output.txt")


## JUST WHITE PHOTOS
## Part 1
## E[y | i %in% I^W] > E[y | i %in% I^B]
##
## where the expectation is over all respondents assigned to M1 or M2,
## all positions in S, and all photos in J^B

mydata <- read.csv("ScaleRaceSpring2017clean.csv")

## just black respondents
mydata.b <- mydata[mydata$R.race == "Black",]
mydata.b <- mydata.b[!is.na(mydata.b$R.race),]
## just white respondents
mydata.w <- mydata[mydata$R.race == "White",]
mydata.w <- mydata.w[!is.na(mydata.w$R.race),]

## keep all photos but nothing else
inds <- 1:24
keep.vars <- c(paste("X", inds, "_Q339", sep=""),
               paste("X", inds, "_Q329", sep=""),
               paste("X", inds, "a.id", sep=""))

mydata.b <- mydata.b[, keep.vars]
mydata.w <- mydata.w[, keep.vars]
## convert factor to character
for (i in 49:72){
    mydata.b[,i] <- as.character(mydata.b[,i])
    mydata.w[,i] <- as.character(mydata.w[,i])    
}
## get rid of black photos
for (i in 1:24){
    i2 <- i+24
    i3 <- i+48
    indic.b <- grepl("mb", mydata.b[,i3])
    indic.w <- grepl("mb", mydata.w[,i3])
    mydata.b[indic.b, i] <- NA
    mydata.w[indic.w, i] <- NA    
    mydata.b[indic.b, i2] <- NA
    mydata.w[indic.w, i2] <- NA    
}
mydata.b <- mydata.b[,1:48]
mydata.w <- mydata.w[,1:48]

## reshape into long format
mydata.b.long <- reshape(mydata.b, direction="long", varying=1:48, v.names="MM")
mydata.w.long <- reshape(mydata.w, direction="long", varying=1:48, v.names="MM")

## keep only obs with non-missing values for MM
mydata.b.long <- na.omit(mydata.b.long)
mydata.w.long <- na.omit(mydata.w.long)


lm.b.out <- lm(MM ~ 1, data=mydata.b.long)
lm.w.out <- lm(MM ~ 1, data=mydata.w.long)

tab.b <- coef_test(lm.b.out, vcov="CR2", cluster=mydata.b.long$id,
                   test="Satterthwaite")

tab.w <- coef_test(lm.w.out, vcov="CR2", cluster=mydata.w.long$id,
                   test="Satterthwaite")
WminusB.diff <- tab.w[1,1] - tab.b[1,1]
WminusB.var <- tab.w[1,2]^2 + tab.b[1,2]^2
WminusB.SE <- sqrt(WminusB.var)
WminusB.z <- WminusB.diff / WminusB.SE
WminusB.pval <- 2*(1-pnorm(abs(WminusB.z)))

cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("                Averaging Over Only White Photos\n")
cat("----------------------------------------------------------------\n")
cat("Part 1:  E[y | i %in% I^W] - E[y | i %in% I^B] > 0\n")
cat("\n")
cat("----------------------------------------------------------------\n")
cat("Avg. MM score among white respondents\n")
print(as.data.frame(tab.w), digits=4)
cat("\nAvg. MM score among black respondents\n")
print(as.data.frame(tab.b), digits=4)
cat("\n\n----------------------------------------------------------------\n")
cat("WminusB.diff     WminusB.SE     WminusB.z    p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.diff, "       ", WminusB.SE, "    ", WminusB.z, "       ",
    WminusB.pval, "\n")
cat("----------------------------------------------------------------\n")
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")











## Part 2
## E[y | i %in% I^W_{M1}] - E[y | i %in% I^B_{M1}] >
##    E[y | i %in% I^W_{M2}] > E[y | i %in% I^B_{M2}]
## where the expectation is over all respondents assigned to M1 or M2,
## all positions in S, and all photos in J^B

mydata <- read.csv("../Data/Cleaned\ Data/ScaleRaceSpring2017clean.csv")


## just black respondents
mydata.b <- mydata[mydata$R.race == "Black",]
mydata.b <- mydata.b[!is.na(mydata.b$R.race),]
## just white respondents
mydata.w <- mydata[mydata$R.race == "White",]
mydata.w <- mydata.w[!is.na(mydata.w$R.race),]


## keep all photos but nothing else
inds <- 1:24
keep.vars.template <- c(paste("X", inds, "_Q329", sep=""),
                        paste("X", inds, "a.id", sep=""))
keep.vars.notemplate <- c(paste("X", inds, "_Q339", sep=""),
                          paste("X", inds, "a.id", sep=""))

mydata.b.m1 <- mydata.b[, keep.vars.template]
mydata.b.m2 <- mydata.b[, keep.vars.notemplate]
mydata.w.m1 <- mydata.w[, keep.vars.template]
mydata.w.m2 <- mydata.w[, keep.vars.notemplate]

## convert factor to character
for (i in 25:48){
    mydata.b.m1[,i] <- as.character(mydata.b.m1[,i])
    mydata.w.m1[,i] <- as.character(mydata.w.m1[,i])    
    mydata.b.m2[,i] <- as.character(mydata.b.m2[,i])
    mydata.w.m2[,i] <- as.character(mydata.w.m2[,i])    
}
## get rid of black photos
for (i in 1:24){
    i2 <- i+24
    indic.b.m1 <- grepl("mb", mydata.b.m1[,i2])
    indic.w.m1 <- grepl("mb", mydata.w.m1[,i2])
    indic.b.m2 <- grepl("mb", mydata.b.m2[,i2])
    indic.w.m2 <- grepl("mb", mydata.w.m2[,i2])
    mydata.b.m1[indic.b.m1, i] <- NA
    mydata.w.m1[indic.w.m1, i] <- NA    
    mydata.b.m2[indic.b.m2, i] <- NA
    mydata.w.m2[indic.w.m2, i] <- NA    
}
mydata.b.m1 <- mydata.b.m1[,1:24]
mydata.w.m1 <- mydata.w.m1[,1:24]
mydata.b.m2 <- mydata.b.m2[,1:24]
mydata.w.m2 <- mydata.w.m2[,1:24]



## reshape into long format
mydata.b.m1.long <- reshape(mydata.b.m1, direction="long",
                            varying=1:24, v.names="MM")
mydata.b.m2.long <- reshape(mydata.b.m2, direction="long",
                            varying=1:24, v.names="MM")
mydata.w.m1.long <- reshape(mydata.w.m1, direction="long",
                            varying=1:24, v.names="MM")
mydata.w.m2.long <- reshape(mydata.w.m2, direction="long",
                            varying=1:24, v.names="MM")

## keep only obs with non-missing values for MM
mydata.b.m1.long <- na.omit(mydata.b.m1.long)
mydata.b.m2.long <- na.omit(mydata.b.m2.long)
mydata.w.m1.long <- na.omit(mydata.w.m1.long)
mydata.w.m2.long <- na.omit(mydata.w.m2.long)

lm.b.m1.out <- lm(MM ~ 1, data=mydata.b.m1.long)
lm.b.m2.out <- lm(MM ~ 1, data=mydata.b.m2.long)
lm.w.m1.out <- lm(MM ~ 1, data=mydata.w.m1.long)
lm.w.m2.out <- lm(MM ~ 1, data=mydata.w.m2.long)

tab.b.m1 <- coef_test(lm.b.m1.out, vcov="CR2", cluster=mydata.b.m1.long$id,
                   test="Satterthwaite")
tab.b.m2 <- coef_test(lm.b.m2.out, vcov="CR2", cluster=mydata.b.m2.long$id,
                   test="Satterthwaite")

tab.w.m1 <- coef_test(lm.w.m1.out, vcov="CR2", cluster=mydata.w.m1.long$id,
                   test="Satterthwaite")
tab.w.m2 <- coef_test(lm.w.m2.out, vcov="CR2", cluster=mydata.w.m2.long$id,
                   test="Satterthwaite")

WminusB.m1.diff <- tab.w.m1[1,1] - tab.b.m1[1,1]
WminusB.m1.var <- tab.w.m1[1,2]^2 + tab.b.m1[1,2]^2
WminusB.m1.SE <- sqrt(WminusB.m1.var)
WminusB.m1.z <- WminusB.m1.diff / WminusB.m1.SE
WminusB.m1.pval <- 2*(1-pnorm(abs(WminusB.m1.z)))

WminusB.m2.diff <- tab.w.m2[1,1] - tab.b.m2[1,1]
WminusB.m2.var <- tab.w.m2[1,2]^2 + tab.b.m2[1,2]^2
WminusB.m2.SE <- sqrt(WminusB.m2.var)
WminusB.m2.z <- WminusB.m2.diff / WminusB.m2.SE
WminusB.m2.pval <- 2*(1-pnorm(abs(WminusB.m2.z)))

WminusB.diff.diff <- WminusB.m1.diff - WminusB.m2.diff
WminusB.diff.diff.var <- WminusB.m1.var + WminusB.m2.var
WminusB.diff.diff.SE <- sqrt(WminusB.diff.diff.var)
WminusB.diff.diff.z <- WminusB.diff.diff / WminusB.diff.diff.SE
WminusB.diff.diff.pval <- 2*(1-pnorm(abs(WminusB.diff.diff.z)))


cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("                Averaging Over Only White Photos\n")
cat("----------------------------------------------------------------\n")

cat("Part 2: (E [y | i in I^W_M1 ]  -  E [y | i in I^B_M1 ]) -\n")
cat("         (E[y | i in I^W_M2 ]  -  E [y | i in I^B_M2 ]) > 0\n") 
cat("\n")
cat("----------------------------------------------------------------\n")
cat("Avg. MM score among white respondents & no template\n")
print(as.data.frame(tab.w.m1), digits=4)
cat("\nAvg. MM score among black respondents & no template\n")
print(as.data.frame(tab.b.m1), digits=4)
cat("\n\n")
cat("Avg. MM score among white respondents & template\n")
print(as.data.frame(tab.w.m2), digits=4)
cat("\nAvg. MM score among black respondents & template\n")
print(as.data.frame(tab.b.m2), digits=4)

cat("\n\n----------------------------------------------------------------\n")
cat("WminusB.m1.diff  WminusB.m1.SE  WminusB.m1.z   p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.m1.diff, "       ", WminusB.m1.SE, "    ", WminusB.m1.z, "       ",
    WminusB.m1.pval, "\n")
cat("----------------------------------------------------------------\n")

cat("\n\n----------------------------------------------------------------\n")
cat("WminusB.m2.diff  WminusB.m2.SE  WminusB.m2.z   p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.m2.diff, "       ", WminusB.m2.SE, "    ", WminusB.m2.z, "       ",
    WminusB.m2.pval, "\n")
cat("----------------------------------------------------------------\n")

cat("\n\n****************************************************************\n")

cat("----------------------------------------------------------------\n")
cat("WminusB.dif.dif  WminusB.dif.dif  WminusB.dif.dif.z   p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.diff.diff, "       ", WminusB.diff.diff.SE, "    ", WminusB.diff.diff.z, "       ",
    WminusB.diff.diff.pval, "\n")
cat("----------------------------------------------------------------\n")

cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")








sink()

